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SUPE5?S0MIC SEPARATED TURBUI^ENT BOUNDARY- 
LAYER OVER A WAVY WALL 

by 

A. Polak and M.J. Werle 
University of Cincinnati 

SUMMARY 

This research is concerned with the development of a pre— 
fiction method for calculating" detailed distributions of surface 
heating rates ^ pressure and skin friction over a wavy wall in 
a two-dimensional supersonic flow. Of particular interest is 
the flow of thick turbulent boundary layers. The surface 
geometry and the flow conditions considered are such that there 
exists a strong interaction between the viscous and inviscid 
flow. First, using the interacting turbulent-boundary layer 
equations, the problem is formulated in physical coordinates 
and "then a reformulation of the governing equations in terms of 
Levy-Lees variables is given. Next, a nximerical scheme for 
solving interacting boundary layer equations is adapted. A 
number of modifications which led to the improvement of the 
numerical algorithm are discussed. Finally, results are presented 
for flow over a train of up to six waves at various flow 
conditions. Limited comparisons with independent experimental 
and analytical results are also given. 
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NOMENCLATURE 


a 

A 

C 

P 

F 

g 

h 

H 

Ki,K2 

z 

I 

L* 

M 

P 

Pr 

Re 


Re 


t 

T 


U, V 


V 


Amplitude . 

Eddy viscosity dancing function. 

• ^ ^ ^ 2 
Skin friction coefficient, x /p u /2. 

^ CO CO ' 

Constant pressure specific heat. 

Normalized longitudinal velocity, F = u/u . 

Normalized total enthalpy, g = H/H . 

€ 

Heat transfer coefficient. 

Nondixnensional total enthalpy, H = H*/u*^. 

CO 

Constants in eddy viscosity models. 

Viscosity parameter, l=pu/p y 

o s 

Mixing length. 

Reference length. 

Mach number. 

Nondimensional static pressure, p = p*/p* u*^ . 

00 CO 

Prandtl number. 

Turbulent Prandtl number. 

Nondimensional turbulent heat flux rate. 


Reynolds number based on reference viscosity, 


Re^ = J^e„uyy;(uf/C*). 


Reynolds number based on free stream viscosity, 

* * * * 

Re = p u L /y 

CO ^ CO C3 ' *^CO 

Time • 

Nondimensional static temperature, T = T*C*/u*^ 

ca 

Nondimensional x^ and X 2 velocity components, 
u = u /u_^, V = V Re^/"^/u^ . 

Transformed v velocity in the boundary-layer. 
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V 





kinc 


£ 

£ 

n 

0 

0c 

8 


T 

0 

r 

y 

? 

:ri,TT2 


Nondimen si on al coordinates (surface or Cartesian ) , . 

* , * 1/2 * * 

Wavelength. 

u^/T . 
e e 

Pressure gradient parameter, (2C/u ) (du /d?) . 

e e 

Ratio of specific heats, C /C . 

p V . . 

Transverse intezrmittency fianction. 

Nondimensional displacement thickness. 

Incompressible displacement thickness. 

Displacement body height. 

Eddy viscosity. 

Eddy viscosity parameter, e = 1 + ^ r . 

Eddy viscosity parameter, £ = 1 + ^ r . 

Transfoimied normal variable. 

Static temperature ratio, 9 = T/T . 

e 

Surface inclination of the body.- 

Surface inclination of the displacement body. 

Nondimensional momentum thickness. 

Longitudinal intermittency function. 

Nondimensional viscosity, y = y*/y * (u*^/C ') 

r =» p * 

Transformed longitudinal variable. 

Functional grouping in inner region eddy viscosity model. 
Nondimensional turbulent shear stress. 

Nondimensional density, p 


* * 
P /p„ 


Subscripts 


frP. 


Conditions evaluated on the displacement body or at the 
outer edge of the boundary layer. 

F'l-at plate value. 


T7T 



X Index for the longitudinal finite difference mesh, 

w Conditions evaluated at the wall. 

“ Conditions evaluated in the upstream freestream. 

Superscripts 

* Denotes dimensional quantities. 
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INTRODUCTION 


This research is concerned with the two-dimensional supersonic 
flow of thick turbulent boundary layers over a train of relatively 
small wave-like protuberances. Interest in this subject arises 
from the need to predict the extent to which an initially flat 
plate boundary layer has been disturbed by a regular corrugation 
the wall surface. The flow conditions and the geometry 
considered here are such that there exists a strong interaction 
between the viscous and inviscid flow. The problem cannot be 
solved without including interaction effects because classical 
boxmdary layer methods would terminate in a separation point 
singularity. 

To handle the present subject by boundary— layer methods/ a 
technique for treatment of the interacting boundary layer 
equations as well as models for turbulence and for the viscous- 
inviscid interaction process must be available. A numerical 
method for addressing closed bubble separation regions was 
developed by Werle and Vatsa [1] . It was applied to a number of 
laminar separated flow problems including flow over a train of 
sine-wave protuberances [2]. This method uses the interacting 
boundary layer equations with a time-like relaxation concept which 
accounts for the boundary- value nature of the problem. This 
approach is adopted in the present study with the inclusion of 
the eddy viscosity model of Cebeci and Smith into the solution 
scheme. The present form of the numerical algorithm includes 
several modifications to that of the earlier work [ 2 , 3J in order 
to accommodate the turbulent nature of the flow, the thick boundary 
layer, and rhe rather dramatic geometry variations of the wavy v/all. 



It was found that the method was capable of handling the 
interacting turbulent flows of present interest. Solutions 
were obtained for flow of thick turbulent boundary layers over 
a train of waves. The results are presented in terms of surface 
pressure, skin friction and heat transfer distributions. The 
predicted trends are compared with available analytical results 
based on small disturbance theory and with experimental data. 



GOVERNING EQUATIONS 


1. Bomdairy Layer Equations in Physical Coordinates 

The suitability of the interacting boundary layer equations 
for describing the relatively strong streainwise variations in 
the boundary layer characteristics due to sudden changes in the 
body geometry has been, at least for the laminar case, verified 
earlier [1, 2] . This approach is used in the present study in 
which Prandtl's classical bomdary layer equations are adopted 
with the only modification that the pressure variation was not 
prescribed but calculated simultaneously from a viscous-inviscid 
interaction model. 

The boTondary layer approximation in two-dimensional viscous 
flow problems implies that the pressure variation is assumed to 
occur only along one coordinate, taken in the general direction 
of the wall shear layer. The degree of this approximation - 
depends on the choice of the coordinate system. While for very 
thin boundary layers over a corrugated wall, or thick boundary 
layers over a relatively flat wall, surface coordinates v/ere 
suitable, (see Ref. 3) for thick boundary layers flowing over 
a small amplitude wavy wall, Cartesian coordinates were found 
to be more appropriate. Accordingly, the governing equations 
will first be written to apply to both the usual surface 

^ 'k k 

coordinates (s , n ) and the Cartesian coordinates Cx , y ) 

* * 

using the notation Cx^, X 2 ) to denote either of these. Non- 
dimensional variables of order one are now defined according 
to the scheme 



* * 


Xt = x,/L 


1/2 * * 
^2 = Re;:/'' x^/L 


(la) 


u = 


* * 
u /u 


„ 1/2 * . * 

V = Re/ V /u , 

r ' 00 ' 


* * *2 

P = P /pn ; 


* * 

P ^ P /P^ / 


* * . *2 


with 


* * * * *2 * 
Re = u L /y (u ^/c ) 

i 03 CO CO ' p' 


(lb) 


(Ic) 


, * * * * * 

ana u , v , p , p and T represent the mean velocities, pressure, 
density and temperature respectively. 

The turbulent boundary layer equations in these variables 

are ; 

Continuity Equation 
d 


9X^ (PU) + 


(pv) = 0 


( 2 ) 


Momentum Equation 


.'• 18 :"% 


•) = P^U 


du 


U J 

e e dx 


e . 3 , 3u . 

^ sTT <“ aZ: ■" "t> 


(3) 


Energy Equation 


/, 3T 9T 


) = - p U 


du 


e e dx 


u + iH_ (,, 9^ 


■1 


9X, 


9x. 


+ T^) 


a /1L_ ^ 1 

8 X 2 Pr 8 X 2 " 


(4) 


where and are the nondimens ional turbulent stress and turbulent 
heat flux respectively. 

The gas is assumed to be air with constant specific heats and 
constant Prandtl number, Pr = 0.72 v;ith the perfect gas law. 
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state Eauation 


V“1 „ 

P = T (5) 

Boundary Conditions 


u(x^. 

X2> 

= 0 

v(x^ , 

Xj) 

=0 at x^ = x^ (x^) 



w 


Xj) 

- T (xj 
w 1 


and (6) 


u(x^, X2) - Ug 

at X2 “ 

T(x^, X2) = 



where (x- ) describes the body surface contour (x« = 0 in 

w ^ 

surface coordinates, for Cartesian coordinates). 

w 

2 . Turbulence Model 

To obtain closure of the system of equations (2.-6) , models for 
the turbulent stress and turbulent heat flux terms are needed. The 
eddy viscosity concept used in conjunction v/ith Prandtl's mix-incr 
length hypothesis for the^wall layer region is the most widely 
used algebraic model for turbulent stress. A well known repre- 
sentation is the two layer eddy viscosity model of Cebeci and Smith 
which has been very successful in modeling turbulence effects for 
flat plate boijndary layers and other attached boundary layers 
with moderate pressure gradients. Less favorable results are 
obtained when using this model for strongly interacting and separated 
flow regions where it appears to fail conceptually. 



In general turbulent quantities like the Reynolds stress 
are governed by transport equations, thus requiring that the 
turbulence history be accounted for. The eddy viscosity concept 
relates the Reynolds stress to only the local mean flow gradient. 
This corresponds to the physical idea that production of turbulence 
at a point due to interaction with the mean flow is cancelled by 
the dissipation due to its self-interaction (this is referred 
to as the "local equilibrium" concept). In other words, the 
eddy viscosity model is the solution to a truncated transport 
equation. In an effort to better align the predictions for 
separated flows with experimental data, previous investigators 
(see Refs. 4 and 5 for examples) have empirically modified the 
equilibrium eddy viscosity model to account for the history effect. 

’frozen', 'relaxation', and other models were devised and 
successfully applied in several of separated flow predictions. 

One of the present authors [6] also used the 'frozen' and 
relaxation' models in the interacting boundary layer equations 
for separated flows with no significant improvements in the 
predicted results over those obtained with the basic eddy-viscosity 
model. An alternate way to try to achieve more satisfactory results 
would be to model the turbulence via the turbulent transport equa- 
tions, and then eventually introduce modifications which will 
account for the extra effects on turbulence occurring in the 
strongly interacting flows. This will, of course, further tax 
the computing times required for these calculations. 

With the a:bove mentioned limitations in mind, the basic form 
of the Cebeci-Smith eddy viscosity model was adapted for the 
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present study where interaction effects tend to reduce longitudinal 
gradients and only small separation regions are encountered. 

Thus we take 

_ 

^ (7a) 

and relate to by the turbulent Prandtl number as 




(7b) 


The turbulent Prandtl number is here taken constant, Pr^ = 0.90. 

The two layer (outer and inner region) Cebeci-Smith model is then 
given as : 

Inner Region 


U/„>i ^ ^ 

y 


* -*2 , * 
3u 


9x; 


(8a) 


where I* = x* [1 - exp (<-xJ/A*) ] 


(8b) 


with 


== 0,40 and 


* * * * 
A , = 26 (y /p ) (y,, 

w 


3u 


3x, 


* ” 1/2 
/ P ) 


(8c) 


w 


where the absolute value of 3u 73 X 2 has been introduced in 
equation (8c) as a modification of the Cebeci-Smith model for 
reverse flows . 

Outer Re a ion 


(e/m)o = 


* * 

p 


Y 6, . 

2 kinc 


(9a) 



where y is the transverse intermittency function 


Y = {1 - erf[5Cx2/x2 “ 0.78) ] }/2 (9b) 

The variable is the value of x„ at which u/u = 0.995, and 

e / e 

*^kinc incompressible displacement thickness. 


3. Boundary Layer Equations in Transformed Variables 

The boundary layer equations given in Section 1 are here 
recast using the Levy* Lees transformation. 

The new independent variables are defined by 

^ / Pe^e^e '^^l ^ / P "^^2 (10a, b) 


The normalized dependent variables are now defined as: 

velocity ratio, 

F = u/u 
' e 

mean static enthalpy ratio 

ORIGINAL PAGE IS 
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or, mean total enthalpy 

g = h/h^ 

With these definitions, equations C2-4) become: 
Continuity Equation 


IX + 2? — + ■p 

+ 2? — + F 


= 0 


Momentum Eauation 


2C F || + V II 
o ^ 3 t) 


B(e-F^) + |_{ee|X) 

dri 3n 


or 


2^ F — + V 

35 3n 


(It f) Mg-F^) t |^(.I||) 


(lla) 

(llb) 


(11c) 


(12) 


(13a) 


(13b) 
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static Temperature Energy Equation 


25 


39 

95 


+ V 


36 

3n 


aie 


(II) 

^3ti^ 



(A 


g 3 9 t 
Pr 3 ti^ 


or 

Total Tentperature Energy Equation 


2C F 


9g 

35 


+ V 



2a 

2+cs 


[A (g - £/Pr)F 


1_ (ii iS) , 

3n Pr 3n 


(14a) 


(14b) 


where z is the viscosity parameter defined by 

Z = pii/p^y^ (15) 

with p given from Sutherland's viscosity law and the turbulent 
parameters, e and e are defined as . 


g = 1 + (g/p) r 


e = 1 + (s/p) 


Pr 


Pr„ 


(16a) 

(16b) 


where r is the streamwise intermittency function; for fully laminar 
flow r - 0 and for fully turbulent flow r = 1, while for the 

r ' 

transitional region its value varies smoothly from zero to one. 

The parameters a and g are obtained from the local inviscid flow as 


a = u^/T 
e e 

a 

% di 


State Equation 


or 


Pg/p = e 

pVp = V <h 

® u • 
e 



(17a) 

(17b) 

(18a) 

(18b) 


Q 



Boundary Conditions 


P(?,n) 

= 0 

V (?, Tl ) 

= 0 

e (?,n) 


g {?/ n) 



at n = 


(19a) 


where n . = 0 for surface coordinates and r> = n (?) for Cartesian 
” w w 

coordinates. Also we have that 


or 


F(?,n) = 1 

S(?fTi) = 1 as 

g(?/Ti) = 1 


(19b) 


The turbulence relations given in Section 2 can be 
in transformed variables as; 

Inner Region 


(s/ji) 


= 


2 2 „2 2 2 
*^e e ^ 1 ^2 ^1 

‘ I 


/2T 


expressed 


(20a) 


where = 1 - exp(-'ir 2 ) 


(20b) 


p u \i I 

^2 = ^ ^ (e>^ 

26Z B M ^ /2^ 


III ) 

'^iw 


1/2 


(20c) 


Outer Region 


P U Y 5, . 

e 5,6 


(20d) 


5, . 
kxnc 


/2T 


n. 


/Re” p u 
r ^e e 


/ 0 (1”F) dn 


(20e) 
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Note that, as far as the form of the governing equations is con- 
cerned, the only difference between the use of surface coordinates 
and Cartesian coordinates -is in the wall boundary condition 
equation {19a) . This can be eliminated using Prandtl*s trans- 
position theorem by writing that 


? = ? (21a) 

n = n ~ n^(S) (21b) 

V = V - P (21c) 

With these transformation equations (i2-14) and (19) yield; 
Continuity Equation 


3V 

3^i 


+ 2? 


3F 

3? 


+ P = 


0 


(22a) 


Momentum Equation 


2C F — + 
3C 
or 

V M = 

3ii 

0 (8- 

P^) + — (£ 
3n 

- 3F. 
£ — ) 

3n 


(22b) 

2C F ~ + 
3C 

V = 

3n 

(1 + 

j) 0(g-F^) 

+ — (£e 

3n 

3n 

(22c) 

Energy Equation 







2g p ii + 

35 

or 

V = 

3n 

aSLz 

(i£, ' 9_ 

3n 3n 

,l€ 30, 
3i^ 

• 

(22d) 

2| F -^ + 
35 

V M = 

3n 

2a 

2+a 

[£. (s-e/Pr) F 

1- 

3ri 3rj 

(A£ l£) 

(22e) 

Boundary Conditions 







P(5/0) = 

^M5,0) 

- 

0 




6(5/0) = 







or g(5/0) - 

s ^ * 






(23 a) 
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= 1 
= 1 

= 1 (23b) 

The interacting boundary layer calculations require an initial 
velocity and temperature profile at some station ahead of the 
effective interaction region (see Figure 1) . This profile was 
obtained here from a noninteracting two dimensional laminar- 
transitional— turbulent boundary layer calculation bv an ordinary 
marching technique using a prescribed pressure distribution. 


and 

6 (?,«) 
or gU,"^) 


4. Inviscid/Viscous Interaction Model 

The interaction of the boundary layer with the isentropic 
supersonic invrscid flow is modeled in the pressure gradient 
parameter 3 by coupling it to the inclination 9^. The edge 
pressure is obtained from the Prandtl-Meyer relation approximated 
here to second order in terms of 8^ as 


P 


e 



T 


/m?-i 

CO 


(M^ -2)^ + Y^r 

eo ' eo 

4(M^ -1)^ 

CO ' 


where 


6^ = tan 


-1 


(d5^/dx) 


^T " ^w ®S 

0g = tan ^(dx^/dXj^) 


(24) 


(25a) 

(25b) 

(25c) 


6 



I (1 


ou \ , 

e e 


(25d) 
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Once is obtained the isentropic relations and Euler’s equation 
are used to obtain g in equation (17b). Thus, the inviscid and 
viscous flows itiust be solved simultaneously since they are 
directly connected through the displacement thickness given in 
equation (25 c) . 
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NUMERIC-AL METHOD OF SOLUTION 


The numerical method used is an implicit finite difference 
scheme written for the similarity form of the governing eguations 
that marches from some initial station along the surface, to the 
terminal point of interest. To account for the boundary value 
nature of the problem,' Werle and Vatsa [1] have added the time 
dependent concept, similar to the one used for the solution of 
elliptic partial differential eguations. This results in 
modification of only the momentum equation (22b) by replacing 
the pressure parameter 0 with 0 defined as 



This method has been successfully applied to laminar separated 
flow problems with various flow configurations including one with 
multiple interacting regions [-2, 3]. The extension of this 
approach to turbulent boundary layers involves, aside from inclusion 
of the eddy viscosity model into the solution scheme, a number 
of modifications (see also Ref. 6) . Specifically, the fQllowing 
steps were taken. 

1. The numerical stability and convergence rate has been 
enhanced by introducing a new differencing in the continuity 
equation, it has only recently been recognized [7, 8] that the 
longitudinal derivatives in the continuity equation provide a oath 
for interacting flows to propagate information upstream. To 
accommodate this numerically requires the use of some sort of a 
forward difference procedure. In the present work vie ..adopt in 
the continuity equation the following forward differencing 
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F. 

1 


- 

,9P, _ ^i+1 

^ 3 ? i " A? 

where the superscript (o) denotes values at the previous time step 
and subscript i refers to the ith s.tation along the length of the 
surface. 

2. The reliability of the present algorithm was enhanced 
by adopting the 'upwind differencing' concept for the longitudinal 
convection effects. In the reversed flow region upwind differencing 
was used in the longitudinal direction for the convective terms 
in order to satisfy the stability requirements. This eliminates 
the so-called 'artificial convection' concept used earlier [2] 
for the laminar case. This modification is significant because 
the velocities in the reversed flow regions are large'r in the 
turbulent case than in the laminar. Thus the convection term in 
the momentum equation is differenced as 


F. 

X 


(i£) 



+ iF^l) (P^ 


F 


i-1 


)/A? + 



(o) 




(27) 


Pj_ is replaced by for forward flow, and by-F^^°^ for reversed 

flow. By replacing the P^ with F^°^ the occurrence of a separation 
point singularity is avoided [1, 6]. Note that the first term 
on the right hand side of equation (27) vanishes for reversed 
flow, and the second term vanishes for the forward flow. The 
same procedure was followed with the term F30/35 in the energy 
equation. 
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Furthermore, the upwind differencing was fovind also helpful 
in the ri direction, in the convective dominated outer regions of 


the thick turbulent boundary layer. It was brought to our 
* 

attention that upwind differencing of the 3F/3r, term in the 

momentum equation might be required to satisfy the convergence 

criteria of the numerical scheme (see also Ref. 9) . In the 

boundary layer near the wall the diffusion term F dominates the 

• • • nn 

convective- like term F and a central difference scheme for F 

n n 

is appropriate. However, in the outer reaches of the boundary 
layer the diffusion term decreases significantly and numerical 
instabilities occur. From a study of- the model 

equation P + aF = 0 it is found that with central differencina 
nn r) 

the criteria jet An j ±2 must be adhered to, to avoid these 

oscillations. Hence the term F was central differenced when 

!“ i i ^ upwind differenced when j a An j >1. 

3. The convergence rate of the time relaxation solution 

method for the thick boundary layers has been fomd much slower 

than for thin boundary layers . The two cases differ largely in 

that for the thick boundary layer the disturbance of the total 

displacement body from the flat plate value is very small. It 

was argued that the numerical truncation error can be of the same 

order as this relative change per one iteration, thus leading to 

very small convergence rate. We introduced therefore a new 

variable in place of the total displacement body 6^. The. 

is defined as '= [6^(C,t) - ^i ^^i^ 

displacement thickness at the initial station and h is a constant 

s 


* R.T. Davis, Personal Communication. 
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of the order of the maximum protx±>erance height. Calculations 
performed with this modification show improvement in the 
convergence rate. 

The accuracy of the calculated solutions depends on the degree 
of precision of the finite-difference approximation and the step 
size. In turbulent boundary layers, large changes occur in the 
velocity profile in the inner layer very near the surface. A 
sufficient number of mesh points are needed near the wall in order 
to get a good resolution in the predictions of wall shear and 
surface heat transfer. At the outer edge of the bounda 3 ry layer 
where the Levy-Lees variable n acquires large values., the changes 
are, on the other hand, very small. This is especially pertinent 
in the case of a thick turbulent bomdary layer disturbed by a 
relatively small protuberance. Thus, for reasons of efficiency 
and accuracy a variable mesh size in the n direction is used in 
solving most turbulent boxondary layers. A mesh growing in size 
from the wall as a geometric progression is used in the present 
algorithm. Blottner [10] has shown that in terms of a transformed 
normal variable N(n) replacing the stretched Levy-Lees variable ri , 
the truncation error is proportional to AN as N 0, or the 
method of calculation is second order accurate. At the jth grid, 
point the physical coordinate is obtained from 

N./AN 1/AN 

n. = n . (K 3 ° - 1)/(K ° - 1) (28) 

-'max 

where K = Arij/ATij_^, Nj = (j-l)AN, AN = 1, and where AN 

is the constant step in the transformed plane. The second order 
accuracy is achieved by varying j and holding AN fixed [10], 
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It was found here that if instead one replaced AN^ by AN, where AN 

is of course varying with while holding K fixed the error 

diminishes v/ith AN as K i.e. much faster than AN^ . Further 

details are given in Appendix A. Figures 2 show the surface heating 

parameter’s dependence on AN and thus provides an accurate error 

estimation procedure. Based on this step size study it was foxind 

that with values of n = 200, = 55, and K = 1.254, a 7% 

itidA inax 

truncation error was incurred in the calculation of wall heat 
transfer. This represents an acceptable compromise betv/een the 
accuracy and the efficiency of calculations. 

The governing equations were linearized and the partial 
derivatives were replaced by finite differences. The eddy viscosity 
term e/u, appearing as a nonlinear term in the governing equations, 
is approximated by .its previous station value. Central differences 
were used to represent partials with respect to n (except as 
noted above where upwind differencing for F^ was required in the 
outer region of the boimdary layer) as well as for the streamwise 
derivatives of the displacement body height, 6^. Upwind dif ferencincr 
was used in the convective terms in the momentum and energy 
equations and forward differencing with respect to c in the 
continuity equation. 

The calculation commences with certain initial conditions and 
then through the time dependent approach [1] the steady state 
solution for a given set of boundary conditions is: sought. In 
the present calculations the initial conditions were set by taking 
the zero time displacement body to correspond to a flat plate 
boundary layer and the surface protuberance to be of zero height. 
Subsequent time sweeps are conducted with the wave amplitude 
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increasing gradually by a small amount. After the desired geometry 
is reached (after the first 10-15 sweeps) the time-like relaxation 
process is continued until the flow properties are relaxed to 
the final value. This process is shown in Figure 3a where 
the skin friction coefficient at one location (s = 3.58) is shown 
as a function of time iteration n\rnber. This location is near 
the junction of the flat plate with a single sine— wave protuberance 
where separation occurs. The resulting skin friction and 
surface heating distributions are shown in Figure 3b. For this 
case with a thin boundary layer,, the calculation was performed 
in surface coordinates. Figure 3a shows that once the full 
protuberance height is attained (11 sweeps) it takes about 50 
more sweeps for the skin friction to attain its ’steady state' 
value. This calculation, with 41 normal grid points and 71 
longitudinal grid points were performed in 5 minutes of 
computer time on the IBM 370-168. 



RESULTS AND DISCUSSION 


A major interest of the present investigation is in the 

numerical predictions for thick turbulent boundary layers over 

a wavy wall, as those in the experiments of reference [11] . 

The geometry and the flow conditions were therefore chosen to 

coincide with those given in reference [11] . The amplitude 

* * 

and wave length are a =0.29 cm, w = 3.66 cm respectively. 

* 

A reference length L = 15.25 cm was chosen. The base flow 

g 

conditions are defined by = 2.53, Re^ = 10.82 x 10 /m, 

T^ = 174 °K and tVt^ = 0.81. Henceforth, v/e refer to these 
» w o 

conditions as standard flow conditions. 


To obtain the present results it is first necessary to 
generate initial profiles at some point ahead of the first 
protuberance-flat plate jimcture. For the standard flow condi- 
tions this station was taken at x = 72.90, where the initial 
profiles were obtained from a noninteracting calculation to 
correspond to the boundary layer as it develops along the wall 
of the UPWr Langley Wind Tunnel [11] . The details concerning 
the calculation of the initial profiles are described in 
Appendix B. The interacting algorithm v/as subsequently employed 
between this initial station and a dov/nstream station past the 
last protuberance. The problem was first formulated and solved 
in the customary surface coordinates. It turned out that the 
geometry extremes make the use of the Cartesian coordinates 
version of the boundary layer equations more reasonable. The 
results of the calculations shown here were performed with a 


longitudinal stepsize Ax 
boundary layer. 


0.02, and a 55 point grid across the 
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Examples from the calculated results are presented for flow 
over a train of up to six waves, for Mach numbers M =2.5 and 

CO 

3.5, for Reynolds numbers Re^ = 10.82 x 10 ^/m and 32.46 x 10 ^/m 

and for wall to total temperature ratios T /T =0.40 and 0.81. 

w o 

Figure 4 shows the contour of a train of six v/aves, the 
displacement body and the viscous and inviscid pressure distri- 
ror the standard flow conditions. The difference in the 
inviscid and viscous pressures dramatically shows the effect and 
need for interaction. The pressure is calculated from an 
approximation to the Prandtl-Meyer relation, accurate to second 
order in flow inclination angle. The inviscid pressure is 
calculated using the local body slope, whereas the viscous 
pressure is obtained by using the slope of the displacement 
body (.= 6^ = y^ + o). The difference in the viscous and inviscid 
pressure is due to the difference in amplitudes of the actual (y^) 
and displacement (5^) body. It is interesting to observe that the 
viscous pressure is almost periodic even though the average 
displacement thickness decreases. Figure 5 shows with the 
distribution of pressure the corresponding distribution of 
surface heat transfer and skin friction at the same base flow 
conditions. The pressure peaks and peaks in heating occur at 
about the same location ahead of the body surface peak. The 
peak in skin friction is shifted in the opposite direction. 

^'Jhile i_he pressure distribution is nearly periodic, the heating 
levels _^i^d the skin friction peaks rise in the .dovmstream 
direction. The rate of rise in peak heating is decreasing very 
slowly. These results are in contradistinction to our similar 
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study [3.] of thin laminar flow over a train' of sine-v/aves , where 
the peaks in heating decreased rapidly in the streamwise direction. 
Figure 6 points out the fact that the local value of the surface 
parameters is almost unaffected by the presence of additional 
downstream disturbance for turbulent boundary layers Ccompars 
also Figures 5 and 6) . Simply, downstream waves have little 
upstream influence and the problem seems localized. Heating 
levels aft of waves grow as the number of waves increases, but 
the downstream skin friction is unaffected by the number of waves « 
To demonstrate the effect on surface properties due to Mach 
number, wall temperature, and Reynolds number, three additional 
cases are shown in Figures 7A,B through 9A,B. The increase of 
Mach number (Figures 7A,B). from 2.5 to 3.5 causes a decrease 
in the ratio of h /h^ . As in the standard flow case, 

the' location which the first wave was placed was chosen in such 
a way that the flat plate boundary layer displacement 


thickness was about the same as in the experimental study of 
reference [11]‘. 

The lowering of v;all temperature to = 0.40 (Figures 8A,B) 

shows a similar trend in h /h^ as for the increase in Mach 

XUaX z • p • 

number. But the absolute rate of surface heating is much higher 

than in the previous case. Interestingly, the h/h^ curve is 

^ * 

smoother here than in other cases. 

Lastly, an increase in Reynolds number, shown in Figures 9A,B 
is seen to cause a slight increase in -the -ratio of peak heating. 


original page is 

^ POOR QUALIOY 


22 



An interesting aspect of the present results is the location 
of the peaks in pressure, heat transfer and shear. The present 
predictions shov7 the peaks in pressure and heating to occur at 
about the same location. This is in agreement with experimental 
observations [11] . The location of the peak pressure in the 
present results is shifted to the right of the location of the 
inviscid peak pressure location (at y 0) by a phase angle of 
about 60°. Theoretical studies by Inger and Williams [12] and 
by Lekoudis et al. [13] predict such a shift. Data from these 
studies given up to = 2.0 show a shift to the left which drops 
off quickly towards zero at = 2. It is therefore possible 
to expect a phase angle in the opposite sense for > 2, as is 
the case in present results . The maximum wall shear location 
obtained from present calculations is shifted to the right of the 
peaks in pressure and surface heating by about 60°. A,ccording 
to theoretical predictions [13] qualitatively such a shift is 
expected. Experimental data available at the same flow conditions 
[11] show a periodic trend in surface pressure as well as in 
the surface heating distribution. The periodic trend in surface 
pressure is observed also in the present predictions. 

Pressure and heating distributions between the second and 
third peak are compared to the experimental data [11], for Mach 
numbers 2.5 and 3.5 in Figure 7B. While the heating dis- 
tributions in the experiments of reference [11] are nearly 
repetditi-ve over consecutive waves, the present predictions 
show a continuous increase over the length of the waves., Note 
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though that in the experimental study there is also indication 
that separation occurs/ while a lack of separation is observed 
in the ^alytical results of Figure 5. The cause of this 
disagreement is not certain but it could well be due to our 
choice of turbulence model or in the fact that the present 
calculations do not simulate well enough all the test conditions 
(three-dimensional effects or boundary layer development on the 
tunnel wall) . Note the calculated boundary layer displacernent 
thickness of the initial profile at station x = 72.90 is 
2.86 cm, close to the value given in reference [11]. However 
the predicted surface heating value at this station is too high 

when compared to experimental data of reference [llj CThe 

. 2 2 
predxcted value is h^^ = 96.5 watts/m ®K. vs 62.5 watts/m °K 

in experiments) . Other prediction methods also typically over 

predict to about the same level the heat transfer rates for 

boundary layers developing along the wind tunnel walls [14] 

thus indicating that some final adjustments may be needed in 

the turbulence model for these flows. 

While the present prediction method has in a certain v;ay 
accounted for the boundary layer displacement effect by using 
the interacting boundary layer equations, the effect of surface 
curvature was neglected. The curvature effect in turbulent 
flows has been summarized in a comprehensive monograph by 
Bradshaw [15] - It has been found that even in cases v/hen the 
longitudinal surface curvature is accounted for -in the governing 
equations the predicted surface aerodynamic parameters Csurface 
heating, wall shear), still deviate considerably from the 
experimental data. Apparently, the streamline curvature has an 



effect on the turbulence structure not included in the existing 
turbulence models. Bradshaw therefore proposed a simple 
correlation scheme applicable to turbulent flows with not too 
large streamline curvature. For the simplest equilibrium model 
used in the present work this correction consists of multiplying 
the mixing length by a factor f = 1 +ae/(3u/Syl, where e = 3v/3x 
is the extra rate of strain induced by the streamline curvature 
(and 3u/3y is the mean rate of strain) . This correction is 
recommended for only 0.5 5 f S 1.5. The constant a is of 
order 10 . Because of the time lag between the first appearance 
of the curvature and its full effect on turbulence the effective 
value of ae is calculated from the lag equation 


106 1^^ = “o® “ 
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(29) 


where 5 is the boundary-” layer thickness, and the constaint is 
taken to be 10. This idea has been implemented here for the 
standard test case by evaluating 3v/3x along a streamline at 
each station x. It was then assumed that this is a representative 
value for all values of y at this station. The 6 was taken to 
correspond to station x = 73.20. The case with base flov/ 
conditions' was recalculated. It was found that a streamline 
had to be chosen in the lower part of the inner layer (passing 
the 21st grid point at x = 73.20) in order not to violate the 
condition 0.5 $ f < 1.5. The results of this calculation are 
shown in .Figures 10A,B. Observe that the peaks in h/h- and 
C- are now lower and slightly dropping along the wall. A 
comparison between the detailed distribution of the surface 
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heating calculated with and without the curvature correction/ 
and against the experimental data is shown in Figure 10b. The 
comparison with e^qjeriments is now more favorable on the 
compression side of the wave. Recall however that in the wind 
t'unnel tests the boundary layer was separated. 

To examine the effect of separation on the character of 
the heating distribution and also to demonstrate the capability 
of the method to solve multiple separation regions, calculations 
without curvature correction were performed with the base 
flow conditions but for a boundary layer developed along a flat 
■plate up to the junction point x = 24.30 (Figures 11A,B) . 

It is noticed (Fig. IIB) that v/hen compared to the unseparated 
case, the heating values are lower on the expansion side of 
the wave. Plow separation occurred in this case due to the 
different history and also, because of the larger pressure 
gradients: the ratio 6/a is here smaller than for the standard 

case and -therefore the slope of the total ‘displacement body 
steeper. This is even more apparent for a, relatively very thin 
boundary layer. Figure 12 shows the results from a calculation 
performed at the same flow condition as in Figure 11, except 
that junction of the flat plate and the first protuberance 
was moved forttfard to x = 3.3 and the wave height was reduced to 
one-half of the previous value. Larger pressure gradients and 
separation regions now appear closer to the center of the 
valleys are observed. Again, the peak heating rates grow in 
the streamwise direction. 
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In reference [11] it was fom:id that unlike for thin bomdary 
layers, the peak heating rates were insensitive to the wave 
amplitude. To examine this question a set of calculations 
were performed for different ratios of wave amplitude to wave 
length, a/w. According to the present predictions the heating 
rates change considerably with the amplitude even for the 
thick boundary layers [Figs. 13A,’B) . For both thick and thin 
boundary layers the peak heating values, (h/hf^p are 

plotted in Figure 14. The present predictions shov 7 similar 
behavior for both thin and thick boundary layers ; a stronger 
than linear increase in peak heating with the increase of 
wave amplitude. Note that the curvature effect on turbulence 
was not accounted for in the results shown in Figures 13 
and 14. The larger curvatures corresponding to. higher amplitude 
waves probably amplify the predicted trend observed in 
Figures 13 and 14. 



CONCLUSIONS 


A n\imerical method capable of handling multiple interactina 
flow regions was adapted to the problem of thick turbulent 
boundary layer over a wavy wall. 

The results of calculations presented in terms of surface 
pressure, skrn friction, and heat transfer distributions disclose 
features distinctly different from the laminar case. The 
present results show a shift in the location of the viscous 
pressure peaks relative to the peaks in the inviscid pressure 
and to the peaks in the wall shear. These phase shifts are 
in qualitative agreement with theoretical predictions based 
on small disturbance theory. • The location of peaks in viscous 
pressure and heat transfer coincide, and the longitudinal 
pressure variation is periodic. This is in agreement with the 
experimental data. The experiments also show periodicity in 
surface heating distribution, while the present results predict 
a continuous increase in heating indicating a possible weakness 
in the turbulence model for a surface with rapidly varying 
curvature. A semi-empirical modification of the eddy viscosity 
model, intended to account for the longitudinal curvature 
effect, aligns the predictions closer to the experimental data. 

It is recommended that future work be carried out with a more 
accurate turbulence model, and. that a more optimal coordinate 
system be adopted for this problem. 
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APPENDIX A 


GRID SIZE EFFECT ON ACCURACY 

The objective is to assess the accuracy of the predicted 
surface bo'undary layer characteristics from the standpoint of the 
normal mesh size. The very thick turbulent boundary layer as 
produced on the wall of the DPTI^T Langley tunnel is of interest 
here. The boundary layer calculation is carried out as an 
non-interacting 2-D laminar- transitional- turbulent boundary 
layer. The corresponding pressure gradient is calculated from 
the Mach number distribution given in Appendix B. 

The free stream conditions in the test section are: M = 

CD 

2.535/ Re /cm = 1.08 x 10^/ T = 314°R, T /T =0.81 and p = 

00 CO ' Q 

199 psf. And Pr = 0.72, Pr^ = 0,9 were taken. The calculation 

commences at s = 0.3 with a laminar boundary layer, which 

becomes fully turbulent at s = 2.1 (As = 0.1). The Mach number 

becomes constant at s = 40.60. The implicit finite-difference 

method (Flugge-Lotz-Blottner-Davis) with non-uiniform grid size, 

varying according to a geometric progression law, is employed. 

The triincation error in this type of calculation is 

proportional to the quantity (ri-,-, - 2n ■ + ri . ) = An . - An ■ , 

J 3”-*- 3 3“1 

and therefore in general not second order accurate. Blottner 
has shown [10] that in terms of a transformed variable, N = N(n) , 
where the coordinate N is obtained by coordinate stretching 
such that AN is uniform, the error is of order AN . Blottner 
substantiates his conclusion by a set of calculations using the 
transformation 

N./AN^ VAN„ 

nj = nj(K ^ ° - 1)/(K ° - 1) , 


where j — 1, 2, J — = number of grid points, 

~ (j“i)/ (J“i) , K — Anj/Art . , , and AN = constant. 

J J o 

A similar type of transformation, corresponding to the 
geometrical variation of An^ is also employed here: 

N./AN 1/AN 

Hj = Hj(K J 1)/(K - 1) 

where AN = 1/(J-1) and K = (= 1.12) does not vary with J. 

Since the local trmcation error is 


E 


j 


‘j + 1 



" 3-1 


^2 

AN^ + 0 (AN^) 
3N 


and 


2 2 
3 n/9N = 


K 


i/AN 


- 1 


‘AN 


InK^) 


N./AN 


K 


we get 


E. ~ (InK ) K 
J o < 


Nj/AN 




1/AN 


- 1 ) 


At a fixed grid location j CN^/AN - j“l) , and for a large 1/AN 
the error should be therefore plotted versus 

A set of calculations was carried out with J = 80, 85, 90, 

95 , 100, 125 and 150 for K = 1.12 and two values of ri _ 

O J 

(rij = 352 and 200) . For the station s = 54.0 (a location on 
the straight section of the tunnel wall) the wall heat 
transfer parameter (3g/3n)^ is shown in Figure 2A. 

error term suggests, the error decreases v/ith decreasing 
^max" plots linearly with as AN ^ 0. It is also 

observed that the two point formula [10] (used in conjunction v;ith 


OF^f‘ is 

WOR QtTALny 


32 



the governing equations evaluated at the wall) for evaluating 
(3g/oTi) gives better results/ especially at the excessive value 

Vax- 

In a second set of calculations the value of K was varied 
together with J in a specified way: 

K = K ^ ; {J-l)m = J - 1 
o o 

Here K = 1.12, m = positive integer, J = 109, n-, - 200. We 

O O u 

take m = 1, 2, 3 or 4 , so that J = 109, 55, 37 or 28. 

Then , 

mN./AN m/AN 


AN = 1/(J~1) = m/(J^“l) = mAN^ 


= 


^ J 


(Jo-l)Nj 


1/AN 

1)/(K^ ° - 1) . 


Thus the transformation is independent of AN and the error 

2 

term, as in reference [10], is now proportional to AN . This 
is exhibited in Fig. 2B where the variation of the surface heat 
transfer parameter at station s = 54.0 is plotted for values 
J = 109/ 55/ 37 and 28 (compare with Figure 2a) . These calculations 
were performed by first setting = 1.12 and J = = 109 

for 0.3 ^s _< 44.3. At s = 44.3, every 2nd, 3rd or 4th grid 
point value was used (corresponding to m = 2 , 3 or 4) and the 
calculation marched (with the new value J = 55 , 37 or 28. and 
K = past station 44.3. 
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In summary, Blottner’s interpretation of the accuracy of 
the variable grid size was extended by showing that the error 
term of the numerical scheme varies linearly with the 
inverse square of the number of grid spacings, (J~l) , when K 
and J are varying in a specified way, but also plots linearly 
with the inverse exponential of (J-1) when K is fixed. 
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APPENDIX B 


INITIAL PROFILES FOR THICK BOUNDARY LAYEP.S 


A major aim of the present research project was to align 
numerical predictions for flow of thick turbulent boundary-layers 
over a wavy wall with the experimental program of reference [llj . 
For this reason a set of calculations were first performed for 
the boimdary layer as it develops along the wall of the UPWT 
Langley Tunnel. The boundary- layer calculation v/as carried out 
as a non-interacting two-dimensional laminar-transitional- 
turbulent boundary- layer developing from' ahead of the nozzle 
throat under a favorable pressure gradient. In the supersonic 
region downstream of the throat the pressure distribution was 
calculated from the sidewall Mach number distribution, using 
isentropic relations. The Mach number distributio n wa s obtained 
from a characteristics net. Prom these data a cubic 
representation for M was assumed of the form 


2 

M = M - K (1 - |-) (1 

-L 



t 


where M^ = test section Mach number at s - s^. The above 

polynomial representation satisfies the condition dM/ds = 0 

at s = s^. At the location of the first characteristic near 

the throat the Mach number was estimated to be M =1.11. The 

a 

measured distance along the wall from s to s^p is 39.6 (- 20 ft). 
At a location s., 8.86 units downstream of s , the Mach num.ber 

P O'-* 

is Mg = 1.84. Letting s^ = 1, Sg = 9.86 and s^ = 40.6. Using 




these values, the constants and can be deteirmined. For 
Ki, K 2 > 0 and K 2 _< 1 (which is the case here) the polynomial 
^®P^®sentation given above yields monotomically increasing 
Mach number. The suosonic-transonic section of the tunnel was 
also assumed to be represented by a cubic polynomial 


M = (1 




with the following properties : 
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At s = s . 

1 

At s = s 

a 


0.3, M 
M = M 


M^ and dM/ds = 0 . 

--i i-(f) 


Three different values for were chosen: 0.01, 0.03 and 0.05. 

It uumed out that the boimdary layer development downstream of 
the throat was not sensitive to these initial values . Takina 
= 0.03, a boundary layer calculation was performed for 
0.3 _< s _< 97.4 with boundary layer profiles punched on IBM cards. 
(At s = 40.6 the Mach number becomes constant on the side wall; 

® the straight section begins.) These profiles were 

used as initial data for the interacting thick turbulent boundary 
layers. The table belov; gives the calculated boundary layer 
displacement thickness distribution along the constant Mach number 
section of the IJPWT Langley timnel wall at ten stations, for 

= 2.53, T^ = 174.4‘>K, T^T^ = 0.81, Re^ = 10.82 x 10^/m, 

Pr = 0.72 and Pr„ = 0.90. 

T 


s 

40.6 

44.6 

48.6 

52.6 

56.6 

60.6 

64.6 

68.6 

72.6 

76,6 

0 (cm) 

1.69 

1.84 

1.98 

2 . 12 

2.26 

2.40 

2.5-4 

2,67 

2.80 

2.93 
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Figure 2 A ACCURACY STUDY OF WALL HEATING LEVEL 
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Figure 3B SURFACE PROPERTIES FOR ONE 
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Figure 4 DISPLACEMENT BODY AND SURFACE PRESSURE - STANDARI 
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Figure 6 DOWNSTREAM EFFECT ON SURFACE 

PROPERTIES 
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Figure 7A MACH NUMBER EFFECT ON SURFACE PROPERTIES 
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Figure 9B REYNOLDS NUMBER 



Figure 10 A CURVATURE CORRECTION EFFECT ON SURFACE PROPERTIE 




Figure 10 B CURVATURE CORRECTION EFFECT ON SURFACE PROPERTIES 







Figure 1 1 B SURFACE PROPERTIES - SEPARATED THICK BOUNDARY LAYE 
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Figure 12 SURFACE PROPERTIES - SEPARATED THIN BOUNDARY LAYER 
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Fiaure 13 A EFFECT OF WAVE HEIGHT ON SURFACE PROPERTIES 
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Figure 13 B EFFECT OF WAVE 





ORIGINAL PAGE IS 
OF POOR QUALITY 



lOa/w 


Figure 14 MAXIMUM SURFACE HEATING AT THIRD PEAK 




